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We introduce and apply a numerically exact method for investigating the real-time dissipative 
dynamics of quantum impurities embedded in a macroscopic environment beyond the weak-coupling 
limit. We focus on the spin-boson Hamiltonian that describes a two-level system interacting with 
a bosonic bath of harmonic oscillators. This model is archetypal for investigating dissipation in 
quantum systems and tunable experimental realizations exist in mesoscopic and cold-atom systems. 
It finds abundant applications in physics ranging from the study of decoherence in quantum com- 
puting and quantum optics to extended dynamical mean-field theory. Starting from the real-time 
Feynman- Vernon path integral, we derive an exact stochastic Schrodinger equation that allows us 
to compute the full spin density matrix and spin-spin correlation functions beyond weak coupling. 
We greatly extend our earlier work (P. P. Orth, A. Imambekov, and K. Le Hur, Phys. Rev. A 
82, 032118 (2010)) by fleshing out the core concepts of the method and by presenting a number of 
interesting applications. Methodologically, we present an analogy between the dissipative dynamics 
of a quantum spin and that of a classical spin in a random magnetic field. This analogy is used to 
recover the well-known non-interacting-blip-approximation in the weak-coupling limit. We explain 
in detail how to compute spin-spin autocorrelation functions. As interesting applications of our 
method, we explore the non-Markovian effects of the initial spin-bath preparation on the dynamics 
of the coherence cr^ {t) and of (t) under a Landau-Zener sweep of the bias field. We also compute 
to a high precision the asymptotic long-time dynamics of cr^(f) without bias and demonstrate the 
wide applicability of our approach by calculating the spin dynamics at non-zero bias and different 
temperatures. 



I. INTRODUCTION 

The coupling of a system to its environment leads to ir- 
reversible energy flovif between system and environment, 
thus giving rise to the phenomenon of dissipation^^''. In 
addition, thermal and quantum fluctuations in the envi- 
ronmental bath cause fluctuations of the system degrees 
of freedom, which results in a Brownian motion^^^. These 
features of bath-induced dissipation and fluctuations oc- 
cur both in classical and quantum systems. A quantum 
system that becomes entangled with the bath exhibits 
decoherence, the suppression of coherence between dif- 
ferent states in the system^'^. The effect of decoherence 
is particularly crucial if one wants to implement a quan- 
tum computer^'^^^'*, where phase coherence between the 
two qubit states |t) and \i) is used as a resource. In fact, 
virtually no system is completely isolated from its sur- 
rounding making dissipation and decoherence ubiquitous 
in physics, chemistry and biology^'^^"^^. 

An impurity spin embedded in a macroscopic environ- 
ment also emerges as an effective model for strongly cor- 
related materials within (extended) dynamical mean-fleld 
theory^^^^'' . In this work, we consider the paradigmatic 
spin-boson modeP'^^, which is a variant of the Caldeira- 
Leggett model^^^^^. Here, the system consists of only 
two states and the environment is described by a bosonic 
bath of harmonic oscillators. 

The spin-boson model with an Ohmic bath is a par- 
ticularly rich model as it exhibits a wealth of interesting 
phenomena such as a delocalization-localization quantum 



phase transition of the spin for sufficiently strong cou- 
pling to the environment^^^'^''. There exist exact map- 
pings to the anisotropic Kondo modeP^^'^*, to the inter- 
acting resonance level model and to the one-dimensional 
Ising model with l/r"^ interaction^^"^" . Various theo- 
retical proposals have been made to experimentally im- 
plement the Ohmic spin-boson model in a controllable 
low-energy circuit. In particular, the recent progress in 
nanotechnology^^""'^ allows for great control on the dis- 
sipation strength of resistive (Ohmic) environments^®^^'* . 
In principle, this development could lead to the real- 
ization of a tunable spin-boson model with microwave 
photons^^"^*^. A tunable spin-boson Hamiltonian with 
Ohmic dissipation can also be realized using trapped 
ions®^ or in cold-atom systems^^"^" , where sound modes 
of a one-dimensional Bose-Einstein condensate mimic the 
bosonic environment. 

Several methods have been devised to investigate the 
dissipative spin dynamics in this system, for exam- 
ple the well-known non-interacting blip approximation 
(NIBA)'''^^ and extensions to it,^^"^'', (non)-Markovian 
master equations, (iterative) path-integral sampling 
and related techniques^^"^^, Monte-Carlo methods on 
the Keldysh contour, ^■'"^^, and various renormalization 
group approaches. ^'''^'^"^"'^^ The Feynman- Vernon real- 
time functional integral formalism^^'^, which will be the 
starting point in the following, is particularly well suited 
to study this class of models since one can easily elimi- 
nate the environmental degrees of freedom. 

Here, our primary goal is to introduce and apply a 
non-perturbative and numerically exact method to in- 
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vestigate the dissipative dynamics of the Ohmic spin- 
boson model beyond the weak spin-bath coupling limit. 
Starting from the real-time functional integral descrip- 
tion, we derive an exact non-perturbative stochastic 
Schrodinger equation (SSE).^^'*"^^^ Compared to earlier 
SSE approaches/^*^"^^"^ our method allows exact consid- 
eration of the initial spin-bath correlations as we derive 
in the main text below. Our approach works both at zero 
and at finite temperatures, and we may easily consider 
a bias field e(t) with arbitrary time-dependence in the 
Hamiltonian. We have previously applied this method 
to investigate the spin dynamics during a Landau-Zener 
sweep e{t) = vt with velocity v.^^^ In addition, it may 
also be applied to other many-body environments, and 
in particular to a fermionic environment, that can be 
represented in the form of a Coulomb gas such as the 
Kondo model-^^'^s-"^ 

The main idea of our method is to recast the prob- 
lem of finding the exact path-integral amplitudes into 
the form of a numerically solvable linear stochastic equa- 
tion^^^'^"'^^'^^®. The quantum spin evolution is given by 
the average solution over different stochastic realizations. 
We explicitly derive an analogy to a classical spin in a 
random magnetic field. Compared to the NIBA, we treat 
the blip-blip interactions in an exact manner here, then 
solve the SSE numerically. 

The starting point to study the influence of the envi- 
ronment on the dynamics of a quantum spin is the spin- 
boson Hamiltonian 

= f + + f E ^''•(^I + bk)+J2 ^'^bW ■ 

k k 

(1) 

The spin is described by Pauli matrices a", a = x,y,z, 
and the operators bk describe bosonic oscillators with mo- 
mentum k and frequency ojfc. They fulfill bosonic com- 
mutation relations [&fc,&g] = Sk.q. In Eq. (1) we have set 
the reduced Planck constant h = 1. The spin part of 
the Hamiltonian contains a tunneling element A, which 
induces transitions between eigenstates {|t) i li)} of tr^. 
It also contains a bias field e{t), which can be time- 
dependent and sets the energy difference between states 
It) and \\.). The spin couples to the bath mode bk via 
the cr^ component and with coupling strength Afc. 

It is well-known that the spin-bath coupling is uniquely 
characterized by the bath spectral function, which we 
assume to be of Ohmic form 

J(w) = TT E ^l^i^ - Wfc) = 27racje~"/"'^ . (2) 

k 

In the following, we take the bath cutoff frequency uic to 
be the largest energy scale in the system. The dimension- 
less parameter a > describes the dissipation strength. 
For a < 1, the interaction with the bath renormalizes the 
characteristic tunneling frequency scale of the spin from 
its bare value A to 

A.^Af-)"^. (3) 



This important energy scale governs the low-energy 
Fermi-liquid fixed point in the delocalized regime a < 1, 
and the spin dynamics for a < l/2'^'^^''^^. At Uc = 1 there 
occurs a localization quantum phase transition where the 
bath completely suppresses tunneling between the two 
spin states, and formally A^ = 0. 

In this article, we are interested to calculate the real- 
time dynamics of the spin (cr"(i)) for different initial 
preparations of spin and bath. We also consider the 
dynamics of the spin-spin correlation function Cz{t) = 
(cr^(i)cr^(O)). We focus on the regime of dissipation 
strength < a < 1/2, where {a^{t)) exhibits damped 
coherent oscillations. We emphasize that perturba- 
tive master-equation approaches fail for a > 0.1 and 
the regime is experimentally accessible'*^'^^"''^. Non- 
perturbative methods are thus required to reliably calcu- 
late the time-evolution of the spin. In this article, we will 
not discuss even stronger spin-bath coupling a > 1/2, 
where the dynamics of {a^{t)) becomes completely in- 
coherent^^^'^^^, before it is completely suppressed for 

Our method allows us to investigate the non- 
Markovian effects of the initial spin-bath preparation on 
the dynamics of (cr"(t)). We distinguish two different 
preparation schemes. In both cases the spin and bath are 
first brought into contact at a time t^. At this time, we 
assume that the spin is in a given pure state, for example 
Ps{to) = It) (tli and the bath is in canonical equihbrium 
at temperature T. The total state at time to is thus 
given by the product state p{tQ) = psito) €5 ps(io) with 
PB{to) - exp(-i/B/T)/Tr[exp(-i/B/T)] where Hb - 
u>kb\bk and the Boltzmann constant is set to = 1. 
We then hold the spin fixed in state |t) until time tj > to. 
This can be achieved, for instance, by applying a large 
bias field e{t) = ea9{-t + tj) with |eo| » A. During the 
time interval [ioi^/] spin and bath are in contact. The 
large-bias constraint is released for t > tj, and the spin 
starts to evolve in time. The system thus starts out from 
a non-equilibrium state at time tj, which is a spin-bath 
product state 

p(ti) = n){t\^PB{ti). (4) 

We now distinguish between two cases: either we send 
Iq — > — oo or we set = tj. This yields different bath 
states pB{ti) at time ti. We study the influence of the 
two different preparation schemes on the spin dynamics 
in detail below. 

We also discuss initial states which are not spin-bath 
product states, and consider a system starting out from 
its equilibrium state at time t — 0. Computing the spin- 
spin correlation function ((T^(i)(T^(0)), we demonstrate 
the effect of initial spin-bath correlations present aXt = 0. 

The structure of the paper is as follows: after this in- 
troduction, we briefly develop the real-time functional 
integral description in Sec. II mainly to introduce our 
notation. In Sec. HI we explain in detail our novel non- 
perturbative stochastic Schrodinger equation approach. 
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We derive all important results that show how to ex- 
actly solve for the driven Ohmic spin-boson dynamics 
provided that Wc ^ A. In Sec. IV, we make an analogy 
of the quantum spin evolution to the dynamics of a clas- 
sical spin in a random magnetic field. We also expose 
the relation to the NIBA in a transparent manner. In 
Sec. V we discuss the influence of the spin-bath prepara- 
tion on the dynamics of the spin. We provide different 
examples where the initial state of the bath has a pro- 
nounced effect on the time evolution of the spin. We show 
that such non-Markovian signatures can be seen in both 
{(J^{t)) and {a^{t)). In Sec. VI, we describe how spin- 
spin correlation functions such as Cz{t) = (cr^(i)cr^(O)) 
can be calculated within the SSE. In Sec. VII, we dis- 
cuss the dynamics of the spin expectation value {a^{t)) 
in various physically relevant situations. This compo- 
nent of the spin is most interesting since, in contrast to 
a^'^{t), it exhibits universal dynamics for large Wc, and 
is thus relevant to the dynamics of a Kondo spin. We 
close this article in Sec. VII with a summary and a dis- 
cussion of open questions and current limitations of the 
SSE method. In the Appendices we provide for complete- 
ness all relevant formulas from the NIBA, corrections to 
the NIBA, and a rigorous Born approximation result of 
Ref. 76. 



II. REAL-TIME FUNCTIONAL DESCRIPTION 

To study the non-equilibrium dynamics of the spin- 
boson model in a non-perturbative way, we employ the 
real-time functional integral description. ■^'^^ In this sec- 
tion, we set the stage and introduce a few technical con- 
cepts following the seminal work of Leggett et al. in 
Ref. 26. Then, in Sec. Ill, in contrast to Ref. 26 we shall 
treat the blip-blip interaction exactly. This affects, for 
example, the long-time behavior of the spin dynamics. 

We are interested to calculate the spin reduced density 
matrix ps{t) — Trspit), where p(t) is the full density 
matrix and Tr^ denotes the trace over the bath. Its 
components can be expressed using real-time functional 
integrals as 

{af\ps{t)\a'f)^ Jva{-) J Va'i-)A[a]A*[a']F[a,a'] , 

(5) 

with crf,o-'f G {|t) j Here, I?cr(-) denotes integration 
over all real-time spin paths a{t) with fixed initial and 
final conditions. A spin path (7{s) jumps back and forth 
between the two values a — ±1. The initial conditions 
describe the preparation of the spin, while different final 
conditions a(t) = aj, (j'(t) = a'^ yield different elements 
of the spin reduced density matrix at time t. 

The integrand in Eq. (5) contains A[a] and yi*[cr'], 
which denote the amplitude of the spin to follow a path 
in the absence of the bath. The effect of the environ- 
ment on the spin is captured by the real-time influence 



functionaP^'^^^ 

F[a,a']^e^p(-- f dsf ds' \-iLi{s ~ s')^{s)7]{s') 
^ ^ Jto J to 

+ L^{s~s')as)i{s')]), (6) 

where we have introduced symmetric and antisymmetric 
spin paths 77(5) — \[<y{s) + <^'{s)\ and £^{s) — \[(j{s) — 
cr'(s)]. It results from an exact integration over the 
bath degrees of freedom"^'^^'^. It contains the real and 
imaginary parts of the force autocorrelation function of 
the environment 'K{X{t)X{Q))T = L2{t) — iLi{t) with 
X^T.k^k{hl + hk) and 

POO 

Li{t) — / duj J {id) sin cut (7) 

^0 

rOO 

L2{t) ^ / duj J (lj) cos ujtcothf3uj/2, (8) 
Jo 

where /3 = l/fesT with temperature T. 

Next, one parametrizes a general (double) spin path 
and inserts it into the functional integral in Eq. (5). Since 
the spin is held fixed at times t <ti, the double spin path 
is constrained to one of the diagonal (or "sojourn" ) states 
{|tt) I lii)}- If we are interested in a diagonal element 
(population) of Ps{t), we fix the final state of the spin 
path to be a "sojourn" state as well. To calculate an off- 
diagonal element (coherence), we let the spin path end 
at time t in an off-diagonal (or "blip") state {\ti} , lit)}- 

For a path that ends in a sojourn state and makes 2n 
transitions at time tj < ti < t2 < ■ ■ ■ < < t along the 
way, we write the spin paths as 

2n 

C(t)=5]s,0(f-t,) (9) 

2n 

v{t)^J2^,e{t~t,). (10) 

The variables {Si,...,S2„} = {6, "Ci; ■ • • : "Cn} with 
= ±1 describe the n off-diagonal or "blip" parts of the 
path spent in the states {|ti) 7 lit)} during times t2m-i < 
t < (m = 1, . . . ,n), where ^(i) = ±1 and r]{t) = 0. 
The variables {Tq, . . . , T2„} = . . . , Vn}, on the 

other hand, characterize the (n -|- 1) diagonal or "so- 
journ" parts of the path during times t2m < t < i2m+i 
(to = 0,...,n), where ri{t) = ±1 and S,{t) = 0. The 
beginning of the initial sojourn is either at tg — >■ —00 
or at to = tj, depending on whether spin and bath 
are in contact at t < tj. We discuss the influence of 
this initial preparation on the dynamics in detail later. 
Formally we have ^2™+! = and the path's bound- 
ary conditions specify 770 and rjn. Altogether, the two- 
spin path is completely characterized by the variables 

{tQ,ti,...,t2n;^l,---,^n]'nO 1,771, ... ,77„„i,77„}. A 

spin path that ends in a blip state is written in an anal- 
ogous way. 
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Using this parametrization of the spin path in Eqs. (9) 
and (10), we may perform the time integrations in the 
influence functional in Eq. (6), which yields 



where 



Qi = exp 



Q2 = exp 



j>k>0 
^ 2n 

~ 51 EjEkQ2{tj - tk) 



j>k>l 



(11) 

(12) 
(13) 



The bath functions Qi.2{t) are the second integrals of 
Li,2{t), i.e. Qi^2 = Li,2- Explicitly, they read for an 
Ohmic spectral density 



Qi{t) ~ 27ratan ^{coct) 
Q2{t) = 7raln(l +a;^t^) +27raln 



\TTt 



(14) 

(15) 



The influence functional is a product of two terms: Qi 
and Q2 ■ While Qi describes a coupling between the blip 
and all previous sojourn parts of the path, the term Q2 
contains the interaction between all blips (including a 
self-interaction) . 

The environment induces a (long-range) interaction be- 
tween the spin path at different times. The state of the 
spin at time t depends on its state at earlier times, which 
leads to a non-Markovian Heisenberg equation of mo- 
tion for the spin. The form of the interaction depends, 
of course, on the spectral density J{uj) and the temper- 
ature T. At zero temperature, for example, one finds 
that L2{t) = 2nauj^{l - ujlt^)/{l + ujlt^f only decays 
algebraically in time. Non-Markovian effects are thus 
pronounced, especially at long times. At high temper- 
atures, on the other hand, the blip-blip interaction be- 
comes short-ranged. In the white- noise limit at T > Wc, 
for example, one derives L2{t) = 2'KakBT5{t) and the 
dynamics is Markovian. 

The path integral of the reduced density matrix in 
Eq. (5) also depends on the free spin-path amplitudes 
A[a] and A* [a']. These amplitudes contribute a factor of 
i^r] A/2 for each transition between a sojourn state 77 and 
a blip state as well as a bias-dependent phase factor 



Hn = exp 



with 



Jti 



(16) 



(17) 



Altogether, the diagonal element of the density matrix 
describing the probability 



to find the system in state |t) at time t is given by a 
series in the tunneling coupling 

"=i '''' 

(19) 

The sum is only over even exponents of A^", because 
we are calculating a diagonal element of ps (t) . The spin 
expectation value {cr^it)) = P{t) can be expressed as 

{a\t))^P{t) = 2p{t)-l. (20) 

In contrast, for an off-diagonal element of ps{t) the path 
ends in a blip state ^2n = ±1 and one finds 



(tiP5(oi;> = ('^+(t)> = *6nE(V) 



iA\2n-l 



X / dt2n-l--- f dh (21) 



ti 



(18) 



where ^2n = 1 for this off-diagonal element and a'^ = 
l(f7^ + ia^). Note the presence of a boundary term in 
Eq. (21) at the final time t, since it now determines the 
end of the last blip, i.e., t2n = t. 

The formal series expansions in Eqs. (19) and (21) are 
exact. What makes these expressions complicated is the 
fact that the coupling between the spin paths in the in- 
fluence functional F is long-range in time. Hence, one 
must consider all terms coupling different blips and so- 
journs. Their analytical evaluation is only possible in 
special cases, e.g. at a = 1/2, or if one simplifies them us- 
ing approximations. The most prominent so-called non- 
interacting blip approximation (NIBA) is discussed in de- 
tail in Appendix A. It simply neglects all interactions be- 
tween different blips. In the next section, we introduce 
a novel method that allows us to take all terms in the 
influence functional exactly into account. In particular, 
unlike the NIBA, we exactly consider the long-range in- 
teractions between different blips. This is achieved by 
mapping the problem onto a linear stochastic equation 
that can be easily solved numerically. 



III. NON-PERTURBATIVE STOCHASTIC 
SCHRODINGER EQUATION METHOD 

We now present a method to evaluate the full spin re- 
duced density matrix ps{t) in a numerically exact man- 
ner. Its element (i| ps{t) \ j) with i,j G {ti4-} is calculated 
by averaging over solutions of a non-perturbative stochas- 
tic Schrodinger equation (SSE). We explicitly derive the 
SSE from the expressions in Eqs. (19) and (21). This 
method works for all temperatures and, importantly, for 
an arbitrary time-dependent bias field e{t). 

In contrast to other numerical approaches such as 
the real-time Monte-Carlo method, ^^''^^'^^^ or the quasi- 
adiabatic path integral scheme, we will not di- 
rectly evaluate the real-time path integral in Eqs. (19) 
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and (21). Instead, we first decouple the terms bilinear in 
the blip and sojourn variables in the influence functional 
Fn = Q1Q2 in Eq- (11) using Hubbard-Stratonovich 
transformations. We then obtain {i\ps{t)\j) as a sta- 
tistical average over solutions of a stochastic Schrodinger 
equation. 

We want to emphasize that our method takes all terms 
in the influence functional exactly into account. In par- 
ticular, we fully account for all interactions between dif- 
ferent blips. Although this method is quite powerful, 
it is so far restricted to the case of an Ohmic bath with 
< a < 1/2 and a large cutoff frequency Wc >• A (scaling 
limit). The reason for this limitation will become clear in 
the following, when we show that Qi greatly simplifies in 
the Ohmic scaling limit. Although we can formally apply 
our approach also for sub-Ohmic and supcr-Ohmic bath 
spectral functions, it remains an open questions to suffi- 
ciently improve the numerical convergence properties to 
make it useful in practice. 



A. Blip-blip interaction part Q2 

We first analyze the Q2 part of the influence functional 
Fn in Eq. (11), that describes the interactions between 
blips. Before we can apply a Hubbard-Stratonovich 
transformation, we must diagonalize the kernel (52 (i) and 
write it in a factorized form as^^^'^^® 



g2(ij -ife) = 7ra Go+ ^ G,„*,„(tj)^'™(tfc) . (22) 

m — 1 



We truncate the sum and keep mmax terms, but always 
check that the final result is independent of mmax- 

To achieve this, we expand Q2(t) in a Fourier series. 
To obtain only negative Fourier coefficients, we rather 
expand Q2iT) = Q2{''') — (32(2) and write 



Q2{t) = Q2(2) 



x/2 



50 



E 

m— 1 



5m COS 



(23) 



on the interval r G (—2,2). Here, t — [t — ti)/ttot is a 
rescaled time which depends on the total length of our 
numerical simulation ttot- The time r = corresponds 
to the initial time t/, when the large-bias constraint on 
the spin is turned-off. The time r = 1 corresponds to the 
final time imax = ^tot + tj of our numerical simulation. 
Note that (32(2) is a constant that depends on the length 
of the simulation ttot- At T = 0, it reads for example 
(^2(2) = 7raln[l -|- Auj^t^^^]. Since we obtain the Fourier 
coefficients gm numerically, this approach is quite general 
and can be used for various forms of the bath correlation 
function Q2{t), as arise for different spectral densities 
J(w). 

From the Fourier expansion, we identify the (trigono- 
metric) eigenfunctions as 5'2fc-i(T) = cos as well as 



*2fe(r) = sin^^. The coefficients in Eq. (22) read 
Go =go + — Q2(2) 

G2/C-I — G2fc — Qk < . 

We can thus write Q2 in factorized form as 
Q2 = exp{-na[ln(l + A^ltl^,) + G] } 



(24) 
(25) 



2n 



(26) 



Z^m=0 



with constant G = 
linim„,x-*oo G = -ln(l -r ■^^c-'tot 
Eq. (26) approaches unity in this limit 



Note that 
4w^t?„t), so the prefactor in 



Eq. (26), we have used that Y^j>k>i'^j^k 



"2 



1. Since G„ 



To derive 

= —n and 

< for TO > 1 it is more appro- 
priate to write ^/aGm = oG^- Next, we decouple 
the blip variables {Sj} in the exponent in Eq. (26) using 
a total of TOmax Hubbard-Stratonovich transformations, 
resulting in 



:p{-na[ln(l+4a.24j+G]} 

/2n 
d5cxp|i^Sj/is(Tj)| . 



(27) 



The integral over the Hubbard-Stratonovich variables 
{sm\ reads 



m—l 



dSr. 



,/2 



/27r 



and we have introduced a real (height) function 



ni—l 



(28) 



(29) 



The function hs (r) contains information about the envi- 
ronment via the eigenfunctions and eigenvalues of the 
bath correlation function Q2{t)- It also depends on 
the Hubbard-Stratonovich variables {sm}, which can 
be interpreted as Gaussian distributed random (noise) 
variables. One thus finds that {hs{t))s = and 
{hs{t)hs{s))s = ctGo — Q2{t — s)/tt, while all higher mo- 
ments vanish. 



B. Blip-sojourn interaction part Qi 

Let us now turn to the Qi part of the influence func- 
tional Fn = Q1Q2, that couples the blip and sojourn 
part of the spin path. It is important to distinguish be- 
tween the flrst sojourn, which occurs during initial time 
to t ^ ti when the spin is polarized, and all other 
sojourns. We thus separate Qi = Qf^'Q^^'^. 
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The contribution of the first sojourn Q^^'^ encodes the 
initial preparation of the system. It is given by the 
terms where k = 0, 1 in Eq. (12) and reads Q^^^ = 
exp[i: Ei=i - to) - Qiitj - h)}]. Within our 

method we can take it into account in an exact manner. 
This is described in Sec. V, and allows us to study the 
effect of the spin-bath preparation on the spin dynam- 
ics. It plays an important role, e.g. for the dynamics of 
(CT^(t)) or if the bias e{t) depends on time. The spin-bath 
preparation is less important for the dynamics of {a^{t)) 
at constant bias""*. 

The contribution of all later sojourns Q^^'' is given by 

the terms with fc > 2 in Eq. (12) and reads Q''^ ^ = 

^^P[l^J2j>k>2^i^kQiitj - tk)]. Fortunately, it takes 
a particularly simple form for an Ohmic bath and if 
A/uic <C 1 and a < 1/2 (scaling limit), where one may 
safely approximate'^^ 



Qi{t) = 27ratan~^(u;ct) « aTr^9{t) . 



(30) 



We thus find that 

2n 

Q^^^ = exp ina ^ S^Tj 

j>k>2 



exp 



n-l 

iira^^k+iVk 

k=l 



(31) 



This reflects the fact that the main contribution to the 
path integral stems from paths with spin-flip separations 
larger than w^T^. 

If we use the scaling form Qi{t) = aTT'^9{t) for the 
first sojourn as well, this corresponds to the spin-bath 
preparation where to = tj. In this case, the complete 
blip-sojourn interaction term Qi = Q^°^ Q^^^ is given by 



Si = exp 



2n 

j>k>Q 



exp 



lira 



^^k+iVk 



k=0 



(32) 



Let us finally note that we can, in principle, deal with the 
blip-sojourn interaction term Qi in a similar way as with 
Q2. In this case, we must first separate the bath corre- 
lation function Qi(t) in the exponent into a symmetric 
part (3i(|t|) and an anti-symmetric part Qi{t), in order 
to extend the sum over the blip and sojourn variables 
to j < k. Then, we can diagonalize the kernels, com- 
plete the square in the exponent and linearize it using 
Hubbard-Stratonovich transformations. The resulting 
expression for the height function /is(t), however, is no 
longer purely real, but also contains an imaginary compo- 
nent. This leads to slow convergence properties, similar 
to the case of the sign problem known from Monte-Carlo 
sampling. This currently limits our SSE approach to an 
Ohmic bath with A/w^ ^ 1 and < a < 1/2 (see also 
Sec. VIII), where Qi{t) = aTT'^e{t). 



C. Stochastic Schrodinger Equation (SSE) 

We now use the form of the influence functional F„ — 
Q1Q2 that we have derived in the last two sections to 
obtain the spin reduced density matrix as a statistical 
average of solutions of a stochastic Schrodinger equation. 
Employing Eqs. (32), (27) and (16) a diagonal entry of 
Psit) can be written as 



p(r) = l+ / dSj2 



n=l 



iAt 



tote 



^{a/2)[ln{l+4uy,^,)+G] N 2« 



dT2n ■ ■■ d.Ti exp 

/O Jo rc 1 

2n 

X ]Jexp[iSj/i(rj)] . 



iira ^ T]k(,k+i 

k=0 



(33) 



(34) 



Here, we have defined the total height function 

h{T) ^ hs{T) + h,{T) . 

It contains both the random height function hsi^r) in 
Eq. (29) as well as the bias-dependent part which reads 



Kir) 



dr'e(T') 



(35) 



Without the summation over the blip and sojourn vari- 
ables {S.jjTjj} the expression in Eq. (33) has the form 
of a time-ordered exponential, averaged over the random 
noise variables {sm}- This summation, however, can eas- 
ily be incorporated into a product of matrices in the vec- 
tor space of two-spin states {|tt) , Iti) , lit) , III)}, that 
read^^'* 



V^Vn 



where 





-i-Ka -ih(T) 






_„-ih(T) 



_gi/i(r) 




gi/l(T) 





_^-iTTa^ih{T) 



(36) 



K, = ^Attot exp{-(a/2) [ln(l + W^t^,) + G] } . (37) 

Note that Vq = \/S.ttot for Wmax — >■ 00. It is worth 
emphasizing that the two-spin basis states simply corre- 
spond to the four elements of the reduced density matrix 
{i\ps \ The final two-spin state withi,j € {t,i} of 
the real-time spin path determines which density matrix 
element {i \ ps{t) \j) is calculated. A product of matrices 
of the type in Eq. (36) automatically satisfies the require- 
ment that transitions between two-spin states occur via 
single spin fiips. A two-spin path consists of an alternat- 
ing sequence of sojourn (diagonal) and blip (off-diagonal) 
parts. The different signs in Eq. (36) stem from the free- 
spin contribution iA^r//2 for each spin flip between blip 
state ^ and sojourn state ry. 
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We finally arrive at the central result of our work. With 
Eq. (36) we can express Eq. (33) as a time-ordered expo- 
nential 

p(t) = J dS{'S?f\Te~'^o'^-^^'^<^>^) . (38) 

Here, T is the usual time-ordering operator. The two- 
spin states and |$/) are the initial and final states 
of the spin path. When calculating the diagonal element 
P{'^) = (tl Psir) It), we thus have <&/ = |tt)- Since we 
consider an initial polarization of the spin in state |t), 
it follows that |$i) = |tt) as well. We can evaluate the 
amplitudes on the right-hand side of Eq. (38) by solving 
the stochastic Schrodinger equation 

|ci,(^)) = v{r) |$(r)) (39) 

with initial and final conditions $ij — (1,0,0,0)"^. The 
vector (1,0,0,0)-^ corresponds to the basis state |tt)- 
The integration J dS over the Hubbard-Stratonovich 
variables is performed by averaging the result over N 
different realizations of the noise variables {s,n}. One 
then obtains p(t) by averaging over the different results 

1 ^ 

^^W-]^E*^''M = ('J'i(^))'S. (40) 

k=l 

where ^i{t) is the first component of \^{t)) and de- 
notes the average over the Hubbard-Stratonovich random 
noise variables {sm}. We note that {^i{t))s is purely 
real as required. The spin expectation value {a^{t)) is 
given by {a'{t)) = P{t) = 2p{t) - 1. 

Other components of the spin reduced density matrix 
PS can be computed simply by using different final con- 
ditions $ J . In order to calculate the off-diagonal element 
(tl Ps{t) It) for instance, we must project onto the final 
state Itt) which corresponds to $/ — (0,1,0,0)-^. In 
this case we also need to consider a boundary term in 
the influence functional at t = T2„ which arises if the 
spin path ends in a blip state"^. It appears as if the 
system steps back to a sojourn state at the final time 
T. This can be implemented by multiplying ^2iT) with 
(A^ftot/(wc/A))", where = A(A/a;c)"/^^"°' is the 
renormalized tunneling element. 

The different spin expectation values are found from 

{a^t)) = 2(A,ttotAK)"($^(ttotT))5 (41) 
(a«(i)) = 2(A,itotAK)"($^'(ttotr))5 (42) 
(a^(t)) =2($'i(itotr))5-l, (43) 

with $^ = Re $Q, = Im and Vq given in Eq. (37). 
Here, we have set tj — for notational clarity. Because 
of the boundary factor (A^ttot A/wc)", the spin expec- 
tation values {a^''^{t)) are not universal functions and 
strictly vanish in the scaling limit A/wc — >■ 0. A universal 
function depends on the bath cut-off frequency Wc only 
through the renormalized tunneling element A^, and is 



thus universal as a function of the dimensionless variable 
y = Art. 

Apart from using the scaling form Qi{t) an'^0{t) 
[see Eq. (30)], the final expressions in Eq. (40)-(43) are 
still exact in the limit Wmax — ^ oo and N ^ oo. In 
practice, of course, we work with finite values of typically 
TOmax ~ 4000 and ~ 10^ - 10^. We always check that 
the final result for p{t) and (t| ps{t) It) is independent of 
Wmax and A'^. As the numerical accuracy of our results 
scales with the number of noise realizations as N^-^^^, 
we are able to routinely calculate the spin expectation 
values ((T^'^^^(t)) up to an (absolute) uncertainty as small 
as 10-*. 



D. Symmetries of the stochastic equations 

It is interesting to note that the differential equations 
in Eq. (39) with the initial condition <I>i = (1,0,0,0)^ in 
fact obey the additional symmetries 

$"(t) = (44) 
$;(t) = $2(r) (45) 
$4(t) = l-$i(r). (46) 

From the eight real variables {<i'i,2,3,47 '^'1,2.3,4} in fact 
only three are independent. Choosing as independent 
variables {^'i, ^'2, ^2} the differential equations (39) read 
explicitly 

^ 2Vo [- sin(/i)$2 + cos(/i)$2] (47) 
dr^2 = Vo[cos{h) sin(7ra) + cos(7rQ;) s'm{h){2^[ - 1)] 

(48) 

dr^'^ = Vq [sin(7rQ!) sm{h) - cos{na) cos{h){2^[ ~ 1)] . 

(49) 

For a particular realization of the random height function 
/i(t), the time evolution described by these equations 
is not unitary. The values of the different components 
{^[,^2,^2} are therefore not bounded. In contrast, we 
can restore the — >■ symmetry by performing two 
simulations with /i(t) = /ie(T) ± hsi^r), which yield the 
two (unbounded) solutions [t). We then take the aver- 
age $(t) = |[«'+(r) + $"(r)], and find that both <^>[{t) 
and ^2{t) are bounded if and only if e = 0. In this 
case, the components are constrained to ^'i{t) G [0,1] 
and ^2{t) e This significantly improves the 

numerical convergence properties, also for the remaining 
component ^2iT) even though it is not bounded. 

IV. ANALOGY TO CLASSICAL SPIN IN 
RANDOM FIELD AND RELATION TO NIBA 

In this section, we show that we can interpret the 
stochastic equation of motion for {a^{t)) = P{t) in the 
scaling limit A/wc ^ 1 and at zero bias e = as that 
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of a classical spin S{t) that rotates in a random mag- 
netic field H{t). Quantum effects and dissipation follow 
from the non-commutativity of rotations around differ- 
ent axes, which are induced by the magnetic field, and 
due to the average over different random magnetic field 
configurations. 

Importantly, however, the quantum-classical analogy 
only applies to the y and z-component of the spin. These 
components are related as {(T^{t)) = — A~^^(cr^(t))'^. 
The analogy cannot be applied to the component {cr^{t))^ 
because the derivation below relies on an assumption that 
is not valid for {a^{t)). In fact, the "classical" Bloch- 
type equation, which we derive in the following, predicts 
(cr^(t)) = which is incorrect. 

We then exploit the analogy to the classical Bloch 
equations to explain the relation of the SSE method 
to the well-known non-interacting blip approximation 
(NIB A). We show that the NIB A follows in our approach 
from neglecting correlations between the classical trajec- 
tory of the spin and the random height function hs{t). 



A. Analogy to classical Bloch equations 



We have previously shown how to obtain the spin 
expectation values {^^'^'''{t)) by solving the stochastic 
equation (39) for different height functions hir) and 
averaging over the results. Different vector compo- 
nents {^a{t))s determine different spin components [see 
Eqs. (41)-(43)]. For a particular realization of the height 
function, however, the components "I'a(i) are unbounded 
and can therefore not be interpreted as components of a 
spin. The situation is different if we explicitly perform 
the sum over the sojourn variables in the expression of 
Qi in Eq. (32), which gives 



Eq. (33) that (cr^(t)) = P{t) = 2p{t) - 1 is given by 



Qi = 



^ (2 cosTra)*^ 



(50) 



For zero bias e = and if we are interested in calculating 
{a^{t)), where the system ends in a sojourn state ?7„, we 
may use that the real-time functional integral expression 
for p{t) in Eq. (19) is invariant under the reversal of the 
sign of all blip variables — >■ — . The imaginary part in 
Eq. (50) does therefore not contribute to {a^{t)), since it 
cancels after summation over ^i. This assumption does 
not hold for {a^{t)), where the system ends in a blip 
state Neglecting the imaginary part in Eq. (50) is 
therefore not justified for {a^{t)). In fact, {cr^{t)) is solely 
determined by the imaginary contribution in Eq. (50), 
which is antisymmetric in ^i. Since we only keep the 
real part in the following, the resulting equations yield 
{a-{t)) ^ 0. 

Keeping only the real part in Eq. (50), we find from 



/oo 



X / dT2n- 

Jo Jo 



0) 
2n 



dTi Y[exp[i^jhs{Tj)] , 
(51) 



with Vq given in Eq. (37). Note that Eq. (5f) only con- 
tains the random part of the height function hs (r) , since 
we have assumed that e = 0. It is worth noting that the 
right-hand side of Eq. (5f)vanishes at the Toulouse point 
a = 1/2. This point marks the boundary between co- 
herent and incoherent dynamics of {a^{t)). For a = 1/2 
the spin-boson model may be solved exactly via mapping 
to a non-interacting resonant level model, and one finds 

(^^(i))a=l/2=exp(-|^t)3^26. 

Since the sojourn variables {?7j} are now absent from 
the expression, we do not need to distinguish between 
the two diagonal states jft) and jH) anymore. It is 
thus sufficient to work with a three-dimensional basis 
{|tt) ) Iti) , lit)}- To perform the summation over the 
blip variables {Sj} via a matrix product, we introduce 
the three-dimensional matrix V3(r). Specifically, we can 
express {a^(t)) as a time-ordered product, averaged over 
the noise variables, as 



(52) 



with the matrix 



V3(r) — VoV2cos 7ra 





gi''s(-!-) 
_g-i/is(r) 



(53) 



The random height function hs{t) is defined in Eq. (29). 
To calculate {a^{t)), we solve the differential equations 



(54) 



with initial and final conditions = (1,0,0)^, and 

take the average over N different noise configurations: 

Similarly to Eq. (39), the differential equations (54) 
with initial condition — (1,0,0)^ obey additional 
symmetries 



Im$i(T) = 

$*(r) = <i>2(r). 



(55) 
(56) 



Only three out of the six variables Re($i,2,3), 
Ini($i 2,3) are thus independent. As before, we choose 
{^[,^'2,^2} = {Re'l'i, Re<l'2, Im$2}, and the differential 
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equations (54) explicitly read 

dr<P'i = 2Vb\/2cos7ra[cos(/is)$2 - sin(/i^)$^] (57) 
dr<^2 = Vb\/2cos7rasin(/is)'I>i (58) 
dr<^>2 = -Vo%/2cos7racos(/i,)$; . (59) 

Denoting a solution for a particular realization of the 
height function /i(t) by $L^''(r), then the solutions for 
[-/i(r)]^are given by <i>'f(T) ^ <i>'+(r), $^-(r) = -$'2+(t) 
and 4>2 (r) = $2^(t). Since we draw the random 
variables {sm} from a Gaussian distribution with mean 
(sm)5 = 0, both /i(r) and [— /i(t)] are equally probable 
and (/i(t))5 = 0. As discussed earlier, Eq. (54) therefore 
(incorrectly) predicts that (cr=^(T)) = 2($2(r))5 = 0. 

A crucial observation is that the system of differential 
equations in Eq. (54) has an integral of motion 

1 = $i(r)2 + 2|$2(t)P - $'i(t)2 + 2[$^(r)2 + $^'(r)2] . 

(60) 

This allows us to introduce an effective classical spin vari- 
able 

S = (5^ 5^ S^-) = {V2<f',,VWl $'0 , (61) 

which is normalized to unit length \S\ — 1. The equa- 
tions of motion for the spin components in Eqs. (57)-(59) 
can be compactly expressed as the classical Block equa- 
tion 

^S{t) = H{t) X S{t) . (62) 

The effective noisy magnetic field H{t) depends on the 
random height function hs{T) and lies in the x-y-plane 

H = HQ{coshs{T),siiihs{T),{)) . (63) 

The amplitude of the magnetic fields reads i/o = 
•\/2VoV2 cos 7ra. The dissipative dynamics of the quan- 
tum spin follows from averaging over different random 
field configurations as 

{a\t)) = {S^'{t))s. (64) 

The quantum problem of the dissipative time evolution 
of the quantum spin component {a^(t)) at zero bias can 
therefore be formulated as an evolution of a classical 
spin S{t) in a random magnetic field H(t). The quan- 
tum nature of the problem is hidden in the fact that 
spin rotations about different axes do not commute and 
through the averaging over different random field con- 
figurations. It is important to keep in mind, however, 
that the quantum-classical correspondence relies on an 
assumption that is not valid for cr^. It is therefore re- 
stricted to the y and z-components of the spin, and pre- 
dicts that {S'={t))s = 0^ (CT^(t)). 



B. Relation between SSE method and NIBA 

We can employ the classical spin description of the pre- 
vious section to make the relation of the SSE method 
to the NIBA transparent. Starting from the "classi- 
cal" Bloch equations of motion in Eq. (62), which yield 
the exact result for (cr^(t)) = {S^{t))s in the limit of 
"m-raaxT N — oo, wc derive an equation that describes 
(cr^(i)) within the NIBA. 

We start from the "classical" Bloch equa- 
tion (62) in the random magnetic field H — 
HQ(^cos[hs{T)],sin[hs{T)],Oy The different spin 

components obey — HyS^, S'^ = —H^S^ and 
= H^S^ — HyS^. For the ^-component we thus 
obtain 

S'{t) = -A^cosina) f ds cos[hsit) ~ hsis)]S' (s) , 
Jo 

(65) 

since Hq = A-y/cos 7ra for rrimax ~^ oo. For a particular 
realization of the noise {sm}, which defines the height 
function hs{t), Eq. (65) is not in convoluted form. In 
general, one cannot write [hs{t) — /is(s)] as a function of 
the time difference (t — s) only. It is thus not possible to 
solve Eq. (65) via Laplace transformation. 

In order to find the time evolution of the quantum spin 
(cr^(i)), we have to average S^(t) over different magnetic 
field configurations 

|(a^(t))^(^^(t))5--A2cos(^a) 

(66) 

It is important to note that there exist correlations be- 
tween the random height function hg (t) and the classical 
spin trajectory S^{s) such that in general 

(^cos[h,{t)-hsis)]S'is))^ 

^ (^cos[hs{t)-h,{s)])^(^S'{s))^. (67) 

These correlations are absent in the initial state at t — 0, 
but are generated over the course of time, as follows from 
the differential equation (65). The correlations are thus 
small at short times t. Note also that since hs{t) ^ y/a, 
the factor cos[hs{t) — hs{s)] ~ 1 for small a ^ 1. At 
a = 0, both classical and quantum spin undergo un- 
damped Rabi oscillations with frequency A. The corre- 
lations between hs(t) and S^(t) thus become more pro- 
nounced for larger values of a. The mean-field decoupling 
anticipated in Eq. (67) can thus be justified at short times 
t and/or small dissipation a. Indeed, we now show that 
one recovers the NIBA from this mean-field decoupling. 
Using the mean- field approximation of Eq. (67), we 
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obtain for the equation of motion of the quantum spin 
d 



dt 



{a^{t)) = {S^{t)), 



-A cos(7ra) 



X j ds {cos\hs{t) — hs{s) 



(68) 



Using the statistical properties of the height function 
{hs{T))s = and {hs{T)hs{s))s = aGo - Q2{t - s)/t:, 
one easily computes the expectation value 

(cos[hs{t) - h,{s)\^ =exp[-g2(i-s)/7r] . (69) 



The equation of motion thus takes the form 
d 



{a'-{t)) = -A^cos(7ra) 



(70) 



X / dsexp[-Q2(^ - s)/Tr] (fT^(s)) 



dt 



which we recognize as the NIBA equation of motion as 
we show in Appendix A 2 following Refs. 128 and 129. 

We note that one also recovers the NIBA equations for 
finite bias eg ^ within the same mean-field decoupling 
scheme of Eq. (67). For non-zero bias, however, one must 
use the equations of motion in Eqs. (47)-(49) to arrive 



26 



at 



-(a^(i))+ dT[f{t-T)+g{t-T){a\T))]=Q, 

(71) 

where f{t-T) = A^e"'?^^*"^)/'' sin(7rQ;) sin[eo(t-T)] and 
git-r) = A2e-'32(*-^)/^cos(7ra)cos[eo(<-T)]. Since the 
norm of the vector S defined in Eq. (61) is not conserved 
for finite ep however, we cannot interpret each trajectory 
$ j [t] for a fixed realization of the noise S as the path of 
a classical spin in a random field. 

To summarize, within our approach we recover the 
NIBA when neglecting the statistical correlations be- 
tween the classical spin trajectories in a random magnetic 
field and the random magnetic field itself (see Eq. (G7)). 
These correlations develop over the course of time and 
thus become more pronounced at longer times. Since the 
magnetic field fluctuations grow with larger dissipation, 
the correlations also increase with a. The derivation of 
the NIBA using our SSE approach thus makes the valid- 
ity of the NIBA at short times and/or weak dissipation 
very transparent. In the following sections, we study the 
differences between these two methods. 



V. SPIN-BATH PREPARATION EFFECTS 

In this section, we investigate the effect of the initial 
spin-bath preparation on the dynamics of the spin. We 
show that the initial preparation can influence the dy- 
namics even at long times. The SSE method is ideally 



suited to consider such situations, because it allows us 
to take the spin-bath preparation exactly into account. 
Specifically, we give two examples: the dynamics of the 
coherence (t) and the behavior of (t) under a linear 
Landau- Zener sweep of the bias. 

We distinguish two different spin-bath preparation 
schemes: in the first one, spin and bath are brought into 
contact at time ~^ — oo, but the spin is held fixed in 
state It) by applying a large bias field until a much later 
time tj, where \ti — to\ ^ oo. By the time i/, the bath 
has relaxed to the shifted canonical equilibrium state 



PB[<yi) = 



exp{-/3[gB + f EfcAfc(6l + &fc)]} 
Tr[exp(-/3[ifB + f Efe + bk)])] ' 



(72) 



with ffi — I corresponding to initial spin state |t). In the 
second preparation scheme, spin and bath are brought 
into contact only at time tg = tj. The initial bath 
state then equals the canonical thermal state, as given 
by Eq. (72) with a^ = 0. 

We first explain how the exact consideration of the 
spin-bath preparation is technically implemented within 
the SSE method, and then investigate two physical situ- 
ations. 



A. Exact consideration of spin-bath preparation 
within SSE method 

The spin-bath preparation is encoded in the contribu- 
tion of the initial sojourn between time to and time ii, 
which denotes the beginning of the first blip. To exactly 
take the spin-bath preparation into account, we must 
therefore use the full form of the bath kernel function 
Qi{t) = 27Tatan~^{ujct) in the contribution of the first 
sojourn. 

The contribution of the first sojourn corresponds to 
the fc = 0, 1 terms in Eq. (12), and explicitly reads 



exp 



^J2^j{Qi{tj - to) - Qiitj - h)} 



(73) 



We can incorporate these terms in an exact manner by 
adding them to the height function /i(t) in Eq. (34), 
which then explicitly depends on ri, the beginning of 
the first blip, and reads 

h{T, Ti) = hs{T) + K{t) - 2a|tan"^ [uScUotiT - n)] 

-5t„,t, tan-i[a;citotT]|. (74) 

The terms hgir) and h^{T) are defined in Eqs. (29) 
and (35). The last term is present only for the prepa- 
ration scheme — ti and absent for tg ^ — oo. 

The fact that the height function /i(t, ti) in Eq. (74) 
now depends on ti forces us to explicitly perform the 
integration over ti g [0, 1] in Eq. (33). We thus randomly 
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pick a uniformly distributed ti G [0, 1] , which determines 
the height function h{T, ri) in Eq. (74). It also determines 
the initial state of the simulation via a single application 
of T/(ti) on = (1,0,0,0)^ as 

1$^,) = -2(0,e'''(^i'^i\-e-'''(^i'^i),0)'^. (75) 

Note that the factors exp(±j7rQ!) that we would ex- 
pect from naively computing — zT^(ri)$i do not occur in 
Eq. (75), because they are taken into account by the last 
two terms in the height function in Eq. (74). 

We then propagate this initial state |$ri) in the inter- 
val [ti, 1] according to Eq. (39), and calculate the proba- 
bility to find the system in state |t) at time t — rttot +ti 
as 

p(r) = l-K(<I>i(T))5. (76) 

From this we find {a^{T)) ~ 2p{t) — 1. The average is 
over N different choices of ti and sets of random variables 
{Sm}- In each individual run, we set |$(r < ri)) = 0, 
since $(t) only accounts for the contribution of paths 
with at least one spin flip. This also explains the differ- 
ence to Eq. (40). The off-diagonal element (t| ps{t) \i) = 
{cr~^) is still determined by ($2)5 and yields the spin ex- 
pectation value {cr^{t)) = 2(ArttotA/a;c)"Re ($2)5 (see 
Eq. (41)). 

B. Dynamics of coherence {a^{t)) 

In this section, we investigate the dynamics of the 
spin expectation value {a^{t)), which describes phase 
coherence between states |t) and The quantum 

(de)coherence can be measured in a persistent current 
experiment or in a SQUID geometry as suggested by 
Refs. 48 and 49. We show that the initial spin-bath 
preparation has a profound influence on its dynamics, 
not just at short times. Of course, in the long-time limit 
the system relaxes to its equilibrium state and {a^{t)) 
becomes independent on the preparation scheme. 

In Fig. 1, we present {a^(t)) for the two preparation 
schemes that we have introduced above. The spin-bath 
interaction is either turned on at to ^ ~ 00 (red line) or 
at to = ti (blue dashed line). While (o''^(t)) clearly ex- 
hibits oscillations if <o ~^ ~oo, it increases monotonously 
for to = tj. The oscillation frequency is of the order of the 
renormalized tunneling element and independent of 
the bath cutoff frequency uJc- For both preparation pro- 
tocols {a^{t)) approaches the correct equilibrium value, 
which we have computed using the thermodynamic Bethe 
ansatz for the interacting resonant level model, which can 
be mapped to the spin-boson model. The spin relax- 
ation to the ground state occurs due to thermalization 
with the bath. 

We can intuitively understand the appearance of the 
oscillations for to ^ in the following way. For this 
spin-bath preparation the initial bath state at t = tj is 
polarized, because the bath has relaxed to the state pb(1) 
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FIG. 1. Coherence (p-^{t)) for different values of a = 
{0.01,0.05,0.1} with UJc = 200A for a = {0.01,0.05} and 
u}c = 50 A for ct — 0.1. Other parameters are A = 1, 
e = 0, mmax = 3000 and iV = 9 x 10^. We present SSE re- 
sults for two different spin-bath preparations: to — > —<x and 
to ~ tj. While {(T^(t)) oscillates for to — >■ —00, it increases 
monotonously for to — tj. We include results of a weak- 
coupling theory beyond NIBA ("NIBA-|-Corr")^'^^'^^, that we 
discuss in Appendix B. SSE curves approach the correct ther- 
modynamic expectation value {a^)oo at long times as calcu- 
lated from Bethe ansatz'^^. We also show {cr'^)oo,Loss/DiVinccnzo 
calculated from a rigorous Born approximation to order a by 
Loss and DiVincenzo^'' . 

in Eq. (72) due to the interaction with the flxed spin. At 
t = t/, all harmonic oscillators arc in the ground state 
of a shifted quadratic potential. This shifted bath state 
acts as a bias fleld e_B(t) in the z-dircction for the spin. 
At t = t/, it reads 

eeiti) = \k{bl + &fe))p,(i) = -2awe . (77) 

k 

Since the spin is released for t > t/, it relaxes towards 
(f^)oo = for e = 0. The polarization of the harmonic 
oscillator bath therefore slowly disappears for t > t/, and 
as a result eB(t) — > 0. For each oscillator the relaxation 
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process occurs on a timescale given by its frequency ujk- 
Due to the presence of many slow modes in the Ohmic 
bath with wj, < A^, the spin thus experiences the bath 
induced bias field eB(t) until times much larger than A"-'^. 

We can understand the oscillations in {a^(t)) by notic- 
ing that the total ( "magnetic" ) field for the spin reads 



B 



(A„0,esW) , 



(78) 



where initially eB{ti) = —2auJc and we have taken into 
account the renormalization of A to A^.. The spin rotates 
around the magnetic field vector B(t) with a frequency 
A sit) = ^y A'^ + esity that is given by the total field 
strength Different components {a°'{t)) are given by 
projections on the different axes. Oscillations in cr^(t) 
thus only occur if the field does not point along the x- 
direction, i.e., only as long as esit) ^ 0. 

It is worth pointing out that although esiti) ^ A^ 
for the parameters in Fig. 1, we observe an oscillation 
frequency of the order of A^, independently of £3(^7), 
which depends on the bath cutoff uic- This can be eas- 
ily understood from the fact that the bath oscillators 
with frequencies uj^. > A^ relax on a fast timescale much 
smaller than A^^. The relevant bath induced detuning 
for t > A~^ is thus rather given by cb ~ —2aAr. 

The amplitude of the oscillations in {a^(t)) is propor- 
tional to the angle 7 between B{t) and the x-axis that 
reads 7 — tanT^ (esit) / A^) . In fact, measuring the os- 
cillation amplitude of {a^{t)) is a new kind of bath spec- 
troscopy. It yields the bath relaxation function esit) in 
Eq. (77), which contains information about the distribu- 
tion of oscillators and their coupling to the spin. 

This clearly non-Markovian effect of the bath initial 
state is captured exactly within SSE, as shown in Fig. 1. 
Since the NIBA yields erroneous results for calculat- 
ing {a^{t)), we compare SSE to predictions of a weak- 
coupling theory beyond the NIBA ( "NIB A-^Corr" ) . This 
approach perturbatively accounts for all interblip cor- 
relations up to first order in a^''^'^-'^^. Details about 
NIBA-|-Corr are provided in Appendix B. 

In the case of to — tj, where the bath initial state 
is unpolarized, SSE and NIBA+Corr agree well at weak 
dissipation (a = 0.01). For slightly stronger dissipation 
of a G {0.05,0.1}, however, the agreement is limited to 
short times only. This reflects the fact that blip-blip 
interactions become more important at longer times. 

Most importantly, in contrast to NIBA-|-Corr, the SSE 
results approach the correct thermodynamic stationary 
value ((T^)oo at long times, independently of the spin-bath 
preparation and for all values of a. The relaxation to the 
ground state occurs since spin and bath thermalize. This 
clearly exemplifies the strength of the SSE approach. In 
Fig. 1, we include the result for {(J^)oo of two different 
calculations. First, using a thermodynamic Bethe ansatz 
for the interacting resonant level modeP^ and, second, 
using a rigorous Born approximation to order a^^ (see 
Appendices B and C for the analytical expressions). 




t (units of A-i) 



FIG. 2. Survival probability p{t) of a free spin for different 
sweep velocities v/A^ — {0.5, 1, 2} and A = 1. After the jump 
at e = 0, p{t) rapidly converges to the classic Landau-Zener 
result piz on a timescale t « A/ v. 



C. Landau-Zener transition 

Another important situation where the spin-bath 
preparation affects the spin dynamics at long times is the 
famous Landau-Zener level crossing problem, where the 
bias field varies linearly in time like e(t) = vt with v > 0. 
Such a Landau-Zener sweep of the bias arises in a variety 
of physical areas such as molecular collisions'^", chemical 
reaction dynamics,'^ molecular nanomagnets,''^' quan- 
tum information and metrology,^'^'''^^^''^'^ and cold-atom 
systems. "'^'^''""'^''^ 

In the absence of dissipation, the Landau-Zener prob- 
lem can be solved exactly. '^°"''*'^ As shown in Fig. 2, the 
(survival) probability p(t) for the spin to remain in its 
initial state |t) shows a jump at the resonance e = at 
t — 0. After the resonance, p{t) quickly converges to its 
asymptotic value 



lim p{t) = exp 

f 00 



7rA2 



Plz 



(79) 



The convergence occurs as soon as e ^ A, i.e., , on a 
timescale t ^ A/v. 

A fundamental question is how the coupling to an en- 
vironment affects the dynamics and the asymptotic value 
of p{t). Analytical results are only know in certain lim- 
j.(-g86, 144-147^ Quite surprisingly, it was proved rigorously 
in Refs. 148 and 149 that at zero temperature, the asymp- 
totic transition probability in the presence of dissipation 
is still given by the classic Landau-Zener result piz , which 
is derived in the absence of dissipation, provided the spin- 
bath coupling is purely longitudinal (i.e., via cr^) and the 
total system is initially prepared in its ground state. The 
proof is valid for any type of bath. 

The proof breaks down, however, for any other ini- 
tial spin-bath state. In general one finds that the initial 
preparation affects the asymptotic long-time value of p{t) 
in the Landau-Zener sweep. Furthermore, the timescale 
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at which the spin reaches its asymptotic hmit is governed 
by the bath cutoff frequency, which is typically orders of 
magnitude larger than A. 

We exemplify this clearly by investigating a model of 
a spin coupled to a single bosonic oscillator mode as 
arises in cavity QED setups^^^"^^'^. It is described by 
the Hamiltonian 



singlc-modc 



= -c7- + '^a' + ^X{b + b^)+uj,b^b. 



(80) 



In Fig. 3, we observe in this toy model that p(t) undergoes 
a sequence of discrete steps separated in time by cuc/v. 
We change the detuning from an initial value of e(t/) = 
— 150 = —6uJc at tj = — 30A^^ to a final value of e(t/) = 
150A = Gwc a.t tf — 30A^^. This behavior occurs since 
the system is driven through a series of avoided crossings, 
which are separated in energy by Wc- Each spin state is 
dressed by a ladder of bosonic states \nb) where ni, G N 
denotes the occupation number of the bosonic mode. 

The probability p{t) only converges towards piz for the 
initial preparation to — > — oo, which corresponds to the 
initial bath state pb(1) a.t t — tj (see Eq. (72) with only 
one bosonic mode here). The timescale of the conver- 
gence is set by the oscillator frequency: t ^ lUc/v, and is 
much larger than A (compare to Fig. 2). If, on the other 
hand, the system does not start out from the ground 
state of the full Hamiltonian in Eq. (80), as is the case for 
t — tj where the initial bath state reads p_b(0) — |0) (0| 
with b |0) = 0, we observe that p{t) does not approach 
Piz at long times. The long-time value of p{t) thus de- 
pends on the initial spin-bath preparation. Of course, 
the difference in the final values depends on the coupling 
strength between spin and bath. 

We have also investigated the Landau-Zener sweep for 
a spin coupled to an Ohmic bath. This is shown in Fig. 4 
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FIG. 3. Toy model result forp(i) under a Landau-Zener sweep 
for V — and different bath initial states. It is only for 
Pfl(l), i.e., to -> ~oo, that p{t) converges towards pi^ at long 
times of 0{uJc/v). Other parameters are lj^ ~ 25A, A = 8A, 
ti = -30A-\ and A = 1. 
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FIG. 4. Survival probability p{t) of spin coupled to an Ohmic 
bath under Landau-Zener sweep. In the upper panel, the 
detuning is swept with velocity v = 5A^ and dissipation 
strength a = 0.2. In the lower panel, velocity is w = lOA^ 
and a — 0.3. Other parameters read ujc = 50A, A = 1, 
mniax = 2000, and ~ lO''. 



for two different velocities v — {5A^,10A^} and spin- 
bath coupling strengths a — {0.2,0.3}. Prominently, we 
find that the jump at resonance at i = is strongly sup- 
pressed due to the coupling to the bath. The size of the 
jump decreases with increasing a. The series of steps, 
which occurred for the single-mode bath, is replaced with 
a smooth decay of p{t) for the continuous bath. The de- 
cay occurs over a timescale governed by the bath cutoff 
frequency Wc- In this decay region, which occurs for in- 
termediate times A/v < t < ujc/v, the spin dynamics is 
universal^ 

At long times, the system converges to the classic 
Landau-Zener result pi^. Interestingly, this is true for 
both initial spin-bath preparations, at least for suffi- 
ciently weak interaction a — 0.2. While the proof 
of Refs. 148 and 149 is only valid for the preparation 
^0 — >■ — oo, we conclude that significant differences in the 
asymptotic limit ofp(t) require large spin-bath couplings, 
at least a > 0.2. This is in agreement with the result for 
a = 0.3 in the lower panel of Fig. 4, where we notice first 
small differences in the long-time value ofp{t). 

In conclusion, the asymptotic long-time value of p{t) 
depends on the spin-bath preparation scheme. Conver- 
gence to the classic Landau-Zener result piz occurs if the 
full system starts out from the ground state, e.g., for 
an infinitely long sweep. Significant differences, however. 



14 



require sufficiently strong spin-bath coupling. 



A. Computation of Sz{t) with SSE method 



VI. CORRELATION FUNCTION 

With the SSE method we can also access the spin-spin 
autocorrelation function 



c.(t) = K(tK(o))T-K) 



(81) 



where cr^(t) is taken in the Heisenberg picture and 
(cr^)T = Tr(cr^e-'''f^)/Tr(e-^^) denotes the equilib- 
rium expectation value at temperature T with respect 
to the full spin-boson Hamiltonian H in Eq. (1). To ob- 
tain limt_>.oo C'z (t) — 0, we subtracted the equilibrium 
value (ct^)t = ^ tanh ^ with A;, = {A'^g + e'^y/^ and 
Aeff = [r(l - 2a)cos(7ra)]i/2(i-")A^. 

In the following, we focus on the symmetric autocorre- 
lation function, which is the real part of Cz{t). The anti- 
symmetric part Xz{t) = —'29{t)l'mC'z{t) can be computed 
within the SSE approach in a similar way. The symmet- 
ric part of the autocorrelation function Sz (t) = Re Cz {t) 
can be expressed as^^^ 



Sz{t)=Ps{t)+ lim Qsit,to), 



(82) 



where Ps{t) is the bias-symmetric part of P{t) = {a^{t)), 
i.e., , Ps{t) — ^[P{e,t) + P{—e,t)]. The remaining part 
Qs{t,to) in Eq. (82) describes the difference between the 
equilibrium autocorrelation function Sz (t) and the single- 
spin expectation value Ps{t) (for zero bias) due to the 
different bath preparation protocols. We always use tj = 
in this section. 

It is worth pointing out that the initial condition for 
P{t = 0) = 1 does not correspond to a small perturba- 
tion. Therefore, P{t) can not be expressed in terms of 
equilibrium correlation functions. The initial spin-bath 
state is a product state p(0) = |t) €5 Pb{^)- In contrast, 
in the case of Sz (t) the system is in its equilibrium state 
at t — 0, where spin-bath correlations are present. These 
equilibrium spin-bath correlations are described by the 
term Q^(i,to) in Eq. (82)i54>i55 

The initial spin-bath correlations lead to significant dif- 
ferences between Sz{t) and P{t), especially at low tem- 
peratures. At T = 0, for example, P{t) decays expo- 
nentially^^^ (or exponentially x power-law^^"') at long 
times (see more later in Sec. VII). In contrast, the long- 
time decay of the autocorrelation is algebraically Sz (t) ^ 
for all a < 1. This includes the exactly solvable 
Toulouse point a — 1/2, where P{t) = exp(— 7i) with 
7 = 7rAV2wc while Sz{t) ~ -{A/Tr-ftf for t oo'^''^^'^. 
We note that the fact that Szit) ~ for a < 1 fol- 
lows very generally from the Shiba relation which 
yields Iiir^^q Sz{uj) = 27ra(xz/2)^|w| ~ a\uj\, where 
Xz = R.exz(w = 0) is the static susceptibility. 



We now show how to calculate Sz{t) using the SSE 
approach. The additional term Qs{t, to) in Eq. (82), that 
describes the spin-bath correlations in the equilibrium 
state reads explicitly^^^ 



^ / ?' A \ 2n+2?n p 

Qs{t,to) = [-tan^ina)] 2^(^j [2cos(7ra) 



n,m— 1 

X / dt2n+2m ' ' ' / 

Jo Jo 

2n+2m 

X - - - 



rtz 
dt2n+l I dt2n ' ' ' / dti 



to 



to 



(83) 



Here, we have assumed a constant bias value e. For sim- 
plicity, we will focus on e = in the following, but it is 
straightforward to include a finite bias into our formal- 
ism. 

If there was no explicit dependence on the first blips at 
negative and positive times {£,i,^n+i}, we could directly 
use the SSE formalism developed for P{t) in Sec. IV to 
calculate Qs{t,to). There we learned that for e = it is 
possible to first perform the sum over sojourn states, and 
define the stochastic Schrodinger equation via a three- 
dimensional matrix V3(t) (see Eq. (53)). In order to keep 
track of the sign of the initial blips at negative and pos- 
itive time, {^i,Cn+i}, we add a fourth state to the SSE 
matrix formalism. For non-zero bias, we would simply in- 
troduce a fifth state in the version of our formalism with 
V4(t). This additional state serves as the initial state of 
the stochastic Schrodinger equation at times t = and 
t = 



|$(i = to)) = |$(i = 0)) = (0,0,0,l)^ 



(84) 



The enlarged matrix is obtained from V^^t) in Eq. (53) 
and now reads 



Vq{t) = 1/oV2cos(to) 

/ e-^''=(^) 








V 







/ 



(85) 



with hs{T) and Vq given in Eqs. (29) and (37). 

We start our simulation at the early time to, which 
should in principle be sent to to — > — oo. In practice, 
we have to take this limit numerically. We use a large 
negative time tg < 0, where \to\ ^ A^^, and check that 
the final result does not depend on to- This is to ensure 
that spin and bath have come to equilibrium by the time 
t = 0. The final time of our simulation is given by imax > 
0. As before, we denote the total length by itot = ^max — 
to, and the rescaled time hy t = {t — io)/itot € [0, 1]. 
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The function Qsit^to) is then calculated from 

Qsit,h) = -tan2(7ra)($i(|to|Atot)>5 

X ($i[(i+|to|)Atot])5, (86) 

where (|io|/itot))5 is the first component of the solu- 
tion at time r = |io|/itot of the system of equations 

^^ Mr)) = Vq Mr)) (87) 

re [OjiolAtot] (88) 
|$(r = 0)) = (0,0,0,1)^, (89) 

averaged over the noise variables S. The function ($i [{t+ 
l^ol)/^tot])5i on the other hand, is the first component of 
the solution at time t = {t + |to|)/itot of the system of 
equations 

^^mr))=VQmr)) (90) 

re [|iolAtot,iAtot] (91) 
|$(r=|io|Atot)) = (0,0,0,1)^, (92) 

again averaged over the noise S. 

Initially, the system is in the sojourn state 
|$(r = 0)) = (0,0,0, 1)"^. Applying Vq once moves the 
system to one of the blip states {|ti) , lit)}- The sign 
of the first blip is taken care of by the different signs 
in the fourth column compared to the first column in 
Eq. (85). After an even number of spin transitions, the 
system has returned to a sojourn state at time t — 
corresponding to r = |to|/itot- Projecting onto the first 
component $i(t = l^ol/^tot) assures that we only con- 
sider paths in Eq. (86) that visit a sojourn state ai t = 0. 
This is required for the symmetric autocorrelation func- 
tion. The sign of the first blip ^n+i at i > is again 
taken care of by starting from state (0, 0, 0, 1)"^ at t = 
and definition of Vq. 

If we were interested in the antisymmetric autocorre- 
lation function Xz{t) = i6'(t)([(T^(t), (T^(0)])t we would 
have to consider paths that visit a blip state at t = 0. 
This quantity has been analyzed within the NRG via a 
mapping to the anisotropic Kondo model in Ref. 158 and 
159. As in the case of the calculation of the coherence 
{a^{t)), if we want to compute Xzit) we are required to 
keep all four two-spin basis states in the SSE calculation 
and use the four-dimensional matrix V4(r). Apart from 
this difference, the calculation of Xzit) is similar to the 
one for Sz{t). 

B. Results for the symmetric autocorrelation 
function Sz 

In Fig. 5, we present results for Sz{t) for different val- 
ues of a. Comparing Sz{t) to P{t) we find quantitative 
differences already within the first oscillation for a > 0.1. 
This is a result of the spin-bath correlations present in 



the initial state, and will become more pronounced at 
longer times. Since we must simulate the dynamics over 
a sufficiently long time interval at negative times [toiO] 
such that the system has reached equilibrium, the com- 
putation of the autocorrelation function Sz{t) using SSE 
is limited to shorter times compared to the computation 
of P(t). 

Still, in Fig. 6 we show the Fourier transform of the 
symmetric autocorrelation function 

Sz{^) = ^ j ^dtSz{t)e'^'. (93) 

for different values of a, where we note that Sz{—t) = 
Sz(t). It exhibits a peak close to the renormalized tun- 
neling frequency A^. The peak is universal, i.e., in- 
dependent of the value of w^. As expected, the peak 
width (height) increases (decreases) with increasing dis- 
sipation strength a. Furthermore, in agreement with the 
Korringa-Shiba relation^'''^"^^^, Sz{io) shows linear be- 
havior Sziui) ~ a\u!\ at low frequencies. 
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FIG. 5. Symmetric autocorrelation function Sz{t) and spin 
expectation value P{t) = {a^{t)} for different values of a = 
{0.1,0.15,0.2,0.25}. Inset shows Qs{t, to) defined in Eq. (83). 
Other parameters read A — 1, tOc = 100 A, e = 0, and 
to — {— 30, — 30, —20, — 20}A~^. we observe a significant 
quantitative difference between Sz{t) and P{t) already over 
the first oscillation period. Within the NIBA both functions 
are predicted to be identical. 
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FIG. 6. Fourier transform of symmetric autocorrelation func- 
tion Sz((^) for different values of a. Other parameters are 
as in Fig. 5. We observe a peak close to the renormalized 
tunneling frequency Ar (independent of uuc- With increasing 
dissipation, the peak width increases while the peak height 
decreases. The low-frequency behavior is Sz{ljj) ~ a\ijj\ as 
expected from the Shiba relation (Refs. 156-159). 



VII. VARIOUS APPLICATIONS 

In this section, we discuss the dynamics of (cr^(i)) in a 
number of different physically relevant situations. It also 
serves to illustrate the capability of the SSE method. 

We first consider the dynamics of (cr^(i)) at zero bias 
e = 0. We confirm the validity of SSE by comparison 
to the NIBA at short to intermediate times and not too 
strong coupling, where the NIBA is valid. We refer to 
Appendix A for details about the NIBA. 

We also closely investigate the long-time behavior of 
(cr^(i)), where the NIBA and corrections to it fail. We 
exploit the fact that we can compute {a^{t)) with great 
numerical precision of 5 x 10~^. Here, we find expo- 
nential decay, possibly with a power-law in the denom- 
inator, in agreement with a non-perturbative prediction 
from conformal field theory^^^ and an expansion around 
the Toulouse point^^^'^^*. 

We then study the dynamics of {a'{t)) for finite static 
bias fields e ^ 0. The SSE method is essentially (numer- 
ically) exact for a < 1/2 in the Ohmic scaling regime 
of A/ujc ^ 1 for any given bias field e(t). This makes 
SSE particularly useful in parameter regions where no 
approximation scheme is known, for instance, at low tem- 
peratures T, small bias fields e and intermediate cou- 
pling strength a > 0.05. In particular, we show that 
a correction to the NIBA for non-zero bias fails already 
for a > 0.05. In Appendix B we include relevant pre- 
dictions of this NIBA correction which we refer to as 
"NIBA+Corr" . 



A. Zero bias dynamics of {a^{t)) 

In this section, we study the dynamics of {a^{t)) for 
zero bias e = 0. We consider both zero and finite tem- 
perature T. We first discuss the dynamics on short- 
to-intermediate timescales, where the NIBA works well. 
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FIG. 7. (a) P{t) = {a^'it)) from SSE for various values of a, 
A = 1, Uc = 100, e = 0, and T = 0. For a given value of 
a, curves corresponding to different values of uoc/ A ^ 1 scale 
on top of each other in units of the renormalized tunneling 
frequency Ar. Numerical parameters are mmax = 3000, A'^ = 
5 X 10^. (b) Quality factor f2/7 of damped oscillations at 
T = 0. Solid line shows prediction from CFT and NIBA. (c- 
d) Finite temperature comparison of {a^{t)) between SSE and 
NIBA at a = 0.1. Other parameters read A = 1, ujc = lOOA, 
e = 0, mmax = 2000, N = 5 x 10"'. We find good agreement 
between SSE and NIBA, quantitative agreement improves for 
higher temperatures as expected. SSE also agrees with NIBA 
temperature T*{a = 0.1) — 2.6A, where the coherent-to- 
incoherent crossover occurs. 



We then consider the long-time dynamics of {a^{t)), 
where NIBA fails. Here, SSE predicts exponential de- 
cay, possibly with a power-law denominator, which is 
in agreement with non-perturbative analytical predic- 
tions^^^'^^'''^^^. The oscillation frequency is found to be 
in excellent agreement with predictions from conformal 
field theory 



1. Dynamics at short-to-mtermediate times 

At zero bias the NIBA is valid for short-to-intcrmediate 
times and not too strong coupling. We refer to Ap- 
pendix A for details. 

At T = 0, the NIBA predicts that (cr'=(t)) is a sum 
of a coherent part Pcoh and an incoherent part -Pine- 
The coherent part describes damped coherent oscillations 
with frequency £7 = Aoff sin 2(1-^0.) ^'^'^ quality factor 
r2/7 = cot 2{i-a) ■ Surprisingly, the same result for the 
quality factor is obtained from a non-perturbative con- 
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formal field theory (CFT) calculation Although the 
NIBA is a weak-coupling approximation, it yields the cor- 
rect quality factor for the full range of < a < 1/2. The 
predicted oscillation frequency, however, is slightly differ- 
ent from the NIBA and CFT. in Fig. 7(a), we present SSE 
results of (o'^(t)) for various values of a. In Fig. 7(b), we 
show that the SSE quality factor precisely matches with 
this formula. 

At finite temperature, the NIBA yields coherent be- 
havior only below a temperature scale T*{oi). The co- 
herent regime is further divided into low temperatures 
T < Aoff and T > Acg. Above T* , the dynamics is 
fully incoherent. In Appendix A 4 we provide all relevant 
NIBA formulas in those parameter regions. In Figs. 7(c) 
and 7(d), we show that SSE agrees well with the NIBA 
over the full temperature range. The quantitative agree- 
ment improves for higher temperatures and for smaller 
values of a (weak-coupling). 



2. Long-time behavior 

Let us now investigate the asymptotic long-time limit 
of (cr^(i)), where |(cr^(i))| < 1. Within NIBA, the alge- 
braically decaying incoherent part Pindt) becomes larger 
than the exponentially decaying coherent part Pcoh(i) af- 
ter a time t » A~g that depends on a. For a = 0.3, for 
instance, one finds that |Pinc| > |-Pcoh| already after one 
half of an oscillation. 

Corrections of the NIBA that take further neighbor 
blip-blip correlations systematically into account modify 
the form of the algebraic power law, but all finite-order 
corrections to the NIBA predict the occurrence of an al- 
gebraically decaying incoherent part. The prediction of 
an algebraic decay of (o'^(t)) at long times is known to be 
an incorrect prediction of the NIBA (and its finite-order 
corrections)'^. 

In contrast to the algebraic decay predicted by the 
NIBA and its corrections, the conformal field theory cal- 
culation in Ref. 125 predicts a purely exponential decay 
of {(J^{t)) at long times. The CFT oscillation frequency 
Als (and decay rate Jls) is also slightly different from 
the NIBA frequency fi, and reads 
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2(1 -a) 
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off 



(94) 



where 

a{a) 



>2(l-a), 



r(i + a)r(l - a) 



l/2(l-a) 



(95) 



In addition, a systematic expansion about the ex- 
actly solvable Toulouse point a — 1/2 — n with 
K <^ 1 yields an exponentially decaying {a^{t)), since 
it yields an incoherent part of the form Pinc(i) = 
-2Kexp[-Aofft/2]/(Aofft)i+2''i24. As shown in Ref. 124, 
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FIG. 8. Long-time behavior of \{a^{t))\ computed with dif- 
ferent methods. Upper panel shows a — 0.25, lower one 
a = 0.3. Other parameters are A = 1, ujc — 50, e = 0, 
T = 0, m^ax = 1000, and AT = 2 X 10^. "NIBA-(-corrections" 
refers to P'^'(Acfft) of Ref. 26, which is a correction to the 
NIBA that takes nearest-neighbor blip-blip correlations into 
account. It improves the NIBA at shorter times, but also fails 
at longer times. "LS" refers to the CFT prediction of Lesage 
and Saleur^^^, and "Toulouse exp." result from an expansion 
around the Toulouse point^^'*. Within numerical accuracy, 
our results are in agreement with exponential decay, possibly 
with a power-law in the denominator^^^'^^*. We can certainly 
exclude purely algebraic contributions. SSE precisely con- 
firms the CFT frequency scale Als- 



a systematic expansion in k shows that interblip corre- 
lations shift the endpoint of the branch cut in (o'^(A)), 
which is responsible for the algebraic decay of in real- 
time within the NIBA, from A = to the non-zero value 
A = — Aoff/2. This behavior is also found in a recent 
study using real-time renormalization grou (RG) and 
functional RG,^^^ where an analytical result for inter- 
mediate times, which is valid to C(|l — 2a|), is reported 
as well. 

In Fig. 8, we present SSE results of |((T^(t))| for a = 
{0.25,0.3} in the long-time limit. We clearly observe 
exponential decay up to a numerical accuracy of about 
5 X 10~^, and no sign of a purely algebraic contribution. 
This agrees with predictions from CFT ("LS") (Ref. 125) 
and the expansion around the Toulouse point ( "Toulouse 
exp.").^^"* It is worth pointing out that SSE oscillations 
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precisely match with the CFT frequency scale in Eq. (94). 
In Fig. 8, we clearly see that the erroneous algebraic term 
dominates the solution of the NIBA and its first-order 
correction ( "NIBA+corrections" ) already after a few os- 
cillations. 



B. Dynamics of {a^{t)) at non-zero bias 

In this section, we discuss the spin dynamics for non- 
zero bias |e| > 0. We focus on the case where |e| ^ 
A. It is well-known that in this case t he NIBA breaks 
down for temperatures below — \/ A^^^ -I- e^s. One 
possibility to go beyond the NIBA is to consider interblip 
correlations up to first order in the spin-bath interaction 
strength a. It is thus limited to weak spin-bath coupling. 
This approach was introduced in Refs. 72 and 73 and 
we provide all relevant results in Appendix A 4 (see e.g. 
Eq. (Bl)). 

In Figs. 9 and 10, we compare SSE to this weak- 
coupling extension of the NIBA ("NIBA-t-Corr") for e — 
—A and two different temperatures T = {10~'^A,A}. We 
find that at low temperatures T = 10"^ A, "NIBA-|-Corr" 
is limited to quite small values of a < 0.01. Even for 
a = 0.05, one observes large differences to SSE. At larger 
temperatures T = A, the damping is much stronger and 
the qualitative agreement improves. This is similar to the 
zero bias case. The overall agreement between SSE and 
"NIBA-|-Corr" at non-zero bias is worse than the agree- 
ment between SSE and NIBA at zero bias. This makes 
the SSE approach a valuable tool to obtain the dynam- 
ics in the presence of non-zero bias, especially at smaller 
temperatures. 



VIII. SUMMARY AND OPEN QUESTIONS 

We want to end with a summary and a discussion of 
a number of open questions related to the SSE method 
and its application to problems beyond the spin-boson 
model. The spin-boson model finds abundant applica- 
tions in physics from quantum computing to the study 
of dissipation induced quantum phase transitions. It is 
realized in a variety of experimental settings most no- 
tably tunable mesoscopic or cold-atom setups. 

In this article, we have exposed in detail a non- 
perturbative numerical method that allows us to exactly 
solve for the spin dynamics {a°'{t)) in the Ohmic spin- 
boson model for a < 1/2. The method can be applied 
provided the bath cutoff frequency is the largest fre- 
quency scale in the problem cjc ^ A. The underlying 
idea of the SSE approach is very general and consists 
of employing Hubbard- Stratonovich identity to trans- 
form the quadratic and time non-local action of the spin 
into a linear and time-local action. The crucial advan- 
tage is that the functional integral over the spin path 
amplitudes can now be exactly calculated by solving a 
linear Schrodinger-type equation. The price to pay is 
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FIG. 9. Comparison of {a^{t)) for non-zero bias e = — A and 
low temperature T = 10~^A between a weak-coupling exten- 
sion of the NIBA ( "NIBA-hCorr" dashed) and SSE (solid). 
Although both solutions agree for a = 0.01, we observe signif- 
icant deviations already for a = 0.05. This is expected since 
the correction to the NIBA takes the interblip correlations 
only up to first-order in a into account. Other parameters 
are coc = 200A, m^ax = 3000 and = 4.5 x 10^. 
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FIG. 10. Comparison of {a^{t)) for non-zero bias e — —A 
and temperature T = A between "NIBA+Corr" (dashed) 
and SSE (solid). For this larger value of temperature, the 
qualitative agreement between the two approaches has im- 
proved compared to Fig. 9. This is expected as the average 
blip length reduces with temperature'^^. Other parameters 
are cuc = 200A, m^ax = 3000 and iV = 4.5 x lO''. 
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the integration over the Gaussian distributed Hubbard- 
Stratonovich variables {s™}. Since the Schrodinger equa- 
tion contains the variables {s,„}, this integration corre- 
sponds to a numerical average over different Schrodinger 
equation solutions. 

The SSE method exhibits very nice convergence prop- 
erties for the Ohmic model with tUc ^ A and a < 1/2, 
since the random height function hs(t) is purely real in 
this case. As a result, each solution of the stochastic 
equation is bounded for zero bias. Even for non-zero 
bias e ^ 0, the individual solutions are well behaved and, 
for example, do not grow exponentially. The situation is 
different, however, for other bath spectral functions such 
as a sub-Ohmic bath. Here, the random height function 
hg{t) acquires an imaginary part, which leads to such bad 
convergence that the approach becomes impracticable. 

We obtain {a°'{t)) as a statistical average over solutions 
of a time-dependent Schrodinger equation, that is easily 
solved numerically by a standard Runge-Kutta solver. 
Therefore, we can easily consider a time-dependent exter- 
nal bias field e{t). Any non-pathological time-dependence 
can be implemented. As an example, in addition to con- 
stant external bias, we have investigated the case of a 
linear Landau-Zener sweep of the detuning. 

Finally, in contrast to earlier stochastic approaches, 
our method allows to take the initial spin-bath prepa- 
ration exactly into account. In particular, a polarized 
bath initial state can have substantial effects on the spin 
dynamics, as we have shown for {a^{t)) for example. 

An interesting further direction is to apply this gen- 
eral idea^^^'^^^ of using Hubbard-Stratonovich transfor- 
mation to obtain a time-local linear action for the impu- 
rity degree of freedom to other impurity problems such 
as the Kondo model, the resonant level model, or the 
Holstein model. The method could also be generalized 
to the case of quantum transport through a quantum 
dot in the Coulomb blockade regime, where quantum 
Monte-Carlo techniques on the real-time Keldysh con- 
tour have been implemented as well.^'^"^^ Another possi- 
ble extension of the SSE formalism is to consider a spin 
with S > 1/2. This increases the number of two-spin 
basis states, that are necessary and this approach is thus 
limited to S* ^ 0(1) in practice. 
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Appendix A: Non-Interacting Blip Approximation 
(NIBA) 

In this section, we derive and discuss the well-known 
Non-Interacting Blip Approximation'''^^. It is essentially 
a short-time and weak-coupling approximation. It be- 
comes exact in the Markovian limit of an Ohmic bath at 
high temperatures. Neglecting blip-blip interactions sim- 
plifies the functional integral expression in such a way 
that it can be solved analytically via Laplace transfor- 
mation. Even if the inverse transformation to real-time 
cannot be performed exactly, much can be learned from 
an investigation of the analytic structure (branch points, 
branch cuts) in Laplace space. 

The NIBA has many shortcomings. Since it is a short- 
time approximation, it always fails at long times (see 
Sec. IV B). In addition, the NIBA can only be used for 
the spin component (cr^(t)). It cannot be used for calcu- 
lating the coherence {(j^(t)) and the spin autocorrelation 
function Sz{t) (except at large temperatures), where it 
fails even at weak spin-bath coupling. For {a^{t)) it gives 
incorrect results at low temperatures and non-zero bias, 
except at very large bias, where is can be justified again. 

In all cases where the NIBA fails blip-blip interactions 
are important. A weak-coupling extension to the NIBA 
that takes blip-blip interactions to first order into account 
in a is derived in Refs. 3, 72, and 73 and is discussed in 
Appendix B. 



1. Derivation from functional integral expression 
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The starting point is the exact expression for the infiu- 
ence functional F„[{Sj}, {Tj}, {t^}] = Q1Q2 with Qi,2 
given in Eqs. (12) and (13) 



Qi = exp 



Q2 = exp 



- ^jTkQi{tj — tk) 

TT ^ — ' 

j>k>0 

^ 2n 

~ X! EjEkQ2{tj - tk) 



]>k>i 



(Al) 
(A2) 



The Qi-part greatly simplifies in the scaling limit 
A/wc for a < 1/2, since one may use 

Qi{t) — Tr^a0(t). Summing over the sojourn variables 
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{t^i, . . . , ?7n-i} results in Eq. (31) 



one finds after rearranging the order of integration 



exp 



ma^^k+iVk =[2cos(7ra)] e . \ \ )/ /^K J 



fc=0 



dU 



(A3) 



For zero bias e = and if we are interested in cal- 
culating {a^(t)), where the system ends in a sojourn 
state rjn, the real-time functional integral expression for 
p{t) in Eq. (19) is invariant under the simultaneous re- 
versal of the sign of all blip variables {^i, ■ ■ ■ , ^n} — >■ 
{— ^1, . . . , — ^„}. Therefore, only the symmetric part 
cos(7ra) of the exponential e"^"^i contributes to {a^{t)), 
and we find Qi = 2"^^ [cos(7rQ;)] . 

The Q2-pa.rt of the influence functional in Eq. (76) 
contains the interactions between all blips. The NIB A 
consists of neglecting all blip-blip interactions apart from 
the blip self-interactions. This relies on the assumption 
that the average time that the system spends in a sojourn 
state is much longer than the average time it spends in 
a blip state^^. The expression of Q2 then becomes 
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(A9) 



where /(A) is the Laplace transform of f{t). The so- 
lution in the time domain is obtained from an inverse 
Laplace transformation via the standard integral along 
the Bromwich contour C^^^ 



2m Jc 



(AlO) 



Even if the inverse transformation cannot be performed 
explicitly, much can be inferred from a study of the ana- 
lytical properties of (ct^(A)), ie., its singularities, branch 
cuts and residua. 



2. Derivation using Heisenberg equations of motion 



A a result, the infiuence functional does not depend on 
the blip and sojourn variables {Cj,'7j} anymore 



Fn[{t,}] = QlQ^ 



NIBA 



2"-i[cos(7ra)]" J|exp 
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(A5) 



The spin dynamics {a^{t)) = P{t) = 2p{t) — 1 follows 
with Hn = 1 for zero bias from Eq. (19) as 

(CT"(t)) = V -A2cos(7ra) / dia™ • • • / dh 
„-n'- Jo Jo 



n=0 
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X Y[ exp 



-Q2{t2j — t2j-l) 



(A6) 



where we have used that } — 2". One can solve 
Eq. (A6) by Laplace transformation^^ 



(a^(A))= / dte-^'{a^{t)). 
Jo 



(A7) 



If we define the function 

f{t) = cos(7ra) exp 



—Q2it) 

TT 



(A8) 



In this section, we present an alternative and physically 
more transparent derivation of the NIBA, which was de- 
rived in Ref. 128. It starts from the polaron transformed 
spin-boson Hamiltonian 

H = U^HU = ^ {a+e'^ + h.c.) + ^ oJkblbk , (All) 

where H is defined in Eq. (1) and the unitary transforma- 
tion reads U = exp{—^a^il) with = —i J^k ^(^fe~^fc)- 
The Heisenberg equation of motion for a^{t) then reads 

CT^(i) = -iAa+(t)e^"(*) -hh.c. (A12) 
It contains ct* (t) which is calculated to 



4W 



^ Jo 



ds aUs)e-"^^'^ 



(A13) 



and crj = (cr+)*. Inserting Eq. (A13) into Eq. (A12) 
yields 

a'{t) = f ds[a"(s)e'"(*)e-'"(") + h.c] . (A14) 

2 Jo 

We now employ two approximations to recover the NIBA. 
First, we assume that the time evolution of the bath op- 
erators is governed by the free bath Hamiltonian Hb = 
'^k^kbj.bk. The reduced density matrix of the bath re- 
mains unperturbed by the spins. Second, we trace out 
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the bath degrees of freedom in a weak-couphng sense by 
writing 

TrB[e^"(*)e-^"(^)] = cxp{^ [*Qi(i - s) ~ Q2{t ^ s)] } , 

(A15) 

which includes the bath correlation functions defined in 
Eqs. (14) and (15). The equation of motion for the spin, 
averaged over the bath, thus becomes 



dsl {a^As)) cos 



Qi(t-s) 



(A16) 



Using the definition of f{t) in Eq. (A8), this can be writ- 
ten as 



(A17) 



(a|(i))+ / ds/(t-.s)(a|(s))=0. 



If we apply a Laplace transformation, we thus recover the 
result for (ct^(A)) within the NIBA, that we have derived 
in the previous Sec. A 1 in Eq. (A9) 



A + /(A) 



(A18) 



3. Zero temperature dynamics 



In this section, we discuss the predictions of the NIBA 
at zero temperature. At T = 0, the Ohmic bath correla- 
tion function reads Q2{t) = 7rQ;ln[l-|-w^i^] . The Laplace 
transform of f{t) in Eq. (A8) is calculated to 



/(A) = Aeff'(Aeff/A)l-2- . 

It contains the effective tunneling element 

r -il/2(l-a) 

Aeff = r(l - 2a) cos(7ra) 



(A19) 



(A20) 



where the renormalized tunneling element is defined as 
A^ = A(A/wc)"/^^""^ • (A21) 

Both Aoff and A^ are smaller than the bare value A, 
because spin transitions are suppressed in the presence 
of a polaronic cloud of bath modes. From Eq. (A9), we 
find that the function ((t^(A)) has a complex conjugate 
pair of simple poles at 



Ai, 



-7 ± ill = Aeff exp 



±i- 



2(1 -a) 



(A22) 



at which point Ai^2 + /(Ai,2) — 0. It also has a branch 
cut along the negative real axis, which ends at the branch 
point A = 0. We note that the branch cut is absent for 
a = and a = 1/2. In the time-domain, the complete 
solution within the NIBA thus reads 



{a'{t)) = P(t) = P,,^{t) + Pi„c(i) 



(A23) 



The poles give rise to damped coherent oscillations of the 
form 



Pcohit) = — ^ — e '^^ cos Vtt 
1 — a 



(A24) 



with frequency ft ~ Aoff sin ^f^-^^^-^ and decay rate 7 — 
Aeff cos 2{\-a) ' quality factor of the oscillations is 

thus independent of Aoff and reads 



— = cot — 

7 2(1 -a) 



(A25) 



The branch cut, on the other hand, yields a (negative) 
incoherent contribution^^ 



„ , , sin27rQ; , 
Pincit)^ / dz 



(A26) 



The incoherent part dominates the dynamics at long 
times Aoffi ^ 1, where it behaves like 



-Pi„c(t) 



1 



(Aoffi) 



2-2a 



This is known to be an incorrect prediction of the NIBA'^. 
Nevertheless, for short to intermediate times, the NIBA 
makes two correct predictions. First that the dynamics 
is universal, i.e., P(t) is a function of the dimensionless 
scaling variable y = Aoffi only. Results of Pit) for differ- 
ent values of Wc collapse on top of each other, if they are 
plotted as a function of y = Aeffi. Second, the quality 
factor of the oscillations in Eq. (A25) exactly agrees with 
results from conformal field theory in the full range of 
< a< 1/2125. 



4. Finite temperature dynamics 

In this section, we discuss the predictions of the NIBA 
at finite temperature. We include this section for com- 
pleteness, it mostly follows Ref. 3. In general, the quality 
of the NIBA improves with increasing temperature, since 
the average blip length decreases for larger temperatures. 
This follows directly from inserting the finite temperature 
bath correlation function 

QJt) = 7raln(l+w^t2) + 27ralnf— sinh— ) (A28) 

into Q^IBA (^4)^ 

The Laplace transform of /(i) at T > is calculated 

^3,26 
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(A29) 



(A30) 



22 



The branch point of /(A)|t=o at A = with the corre- 
sponding branch cut along the negative real axis, that 
occurred at T = 0, turns into an infinite number of sim- 
ple poles at finite temperatures. The poles A„ are located 
on the negative real axis in intervals —^{n + a) < A„ < 

^{—n + a) with n = 1, 2, . . .. The spacing of the poles 
grows linearly with temperature. One can show^'^^, that 
the resulting contribution Pinc(i) = J2n syrp{Xnt) is 
still negative, like the contribution of the branch cut at 
T = 0, and that this part may be neglected for weak 
dissipation. 

The dynamics is dominated by the behavior of the two 
simple poles of {a^{X)) as a function of temperature. The 
poles Xi,2{T) — — 7(T)±iil(T) move toward the negative 
real axis for increasing temperature. At a temperature 
of T = r* (a) they hit the real axis, where 



ioff 



r(a) 



2tt \aT{l ~ a) 



1 + 7racot(7rQ;) + 2\JW{a) 



For weak dissipation a ^ 1, one finds T* 
used the definitions 

M^(q;) — 7racot(7rQ!) — a^.g2(a) 



1/2(1-q) 



(A31) 



We have 



(A32) 



52(«) = 2 [^^'(1 - - ^'(^ + - (A33) 



(7i(a) — a ^ — 7rcot(7ra) . 



(A34) 



Here, "ii^z) is the first derivative of the digamma func- 
tion ■(/'(z)^^'^. The temperature T* separates the regime, 
where Pit) exhibits coherent oscillations {T < T*) from 
the regime, where it exhibits incoherent decay (T > T*). 

Specifically, as long as T < A^g, one finds that 
Ai,2(T) = -7(T) ± in{T) with3 

n{T) = Aeff{l + a[Re V(«/3Aeff/27r) - ln(/3AefF/27r)] } 

(A35) 

7(T) = |aAeff coth(;9Aeff/2) . (A36) 

For larger temperatures T > Aog, but potentially still 
T < T* , one can expand h{\) in Eq. (A29) up to second 
order in A. From A + /(A) = 0, we find the poles Xi^2{T) 
by solving"^ 



(1 — g2u)x^ + {a + giu)x + u = . 



(A37) 



Here, we have introduced x = /3A/27r, u = 
(/3Aoff/27r)^~^"ft,(0). The two poles are complex conju- 
gates Ai.2 = — 7 ± iSl for T < T* . On the other hand, 
for T > T* they are real and negative Ai_2 = ~7i,2 
and lie in the interval (— 27rTa,0). As temperature is 
increased, both poles move in opposite directions. While 
one of them converges for large temperatures toward 



Ai(T Aoff) 0, the other one converges toward 
A2(r> Aefi) ^ -2-KTa. 

The real-time dynamics of the spin is again obtained 
from Laplace inversion. In the coherent regime Agff < 
T < T*, we find that (cr^(t)) exhibits damped oscilla- 
tions with {n, 7} defined by the location of the pole, and 
(f) = — I tan~^(7/r2)| being an initial phase shift. In the 
incoherent regime T > T*, we obtain the dynamics 
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where 71 > 72. Fairly above T* , the prefactor of the 
larger decay rate term, which reads 72/(71—72) becomes 
very small. The decay is thus dominated by the smaller 
decay rate 72, which shows an asymptotic temperature 
dependence of'^ 



72 (T) 



A^/7rT\2c 



2r(a-|-l/2) T VA^y 



2a-l 



(A39) 



In the main text, we compare SSE to these NIBA predic- 
tions and find good agreement at short to intermediate 
times. As expected, the quantitative agreement enhances 
for increasing temperatures. The NIBA fails at longer 
times, non-zero bias fields and off-diagonal elements of 
the reduced density matrix, i.e., {a^{t)). 



Appendix B: Weak-coupling extension to the NIBA 
(NIBA+Corr) 

The NIBA breaks down for finite bias fields e 7^ at 
temperatures below Ah ~ \J A^jj -I- e^^. Even for zero 
bias the NIBA cannot be used for calculating the coher- 
ence {a^{t)). First-order blip-blip interactions are crucial 
in those situations. 

In this section, we state for completeness results of an 
approach that goes beyond the NIBA which was intro- 
duced in Refs. 72 and 73. This weak-coupling extension 
of the NIBA ("NIBA+Corr") considers the interblip cor- 
relations up to first order in the spin-bath interaction 
strength a. It is therefore limited to weak spin-bath cou- 
pling, i.e., small values of a <C 1. 

The first-order interblip contribution to the influence 
functional can be exactly calculated and yields 
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7A2fj 7^Po< 
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with non-zero long time value 
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and 
A, 



= a/A 



(B3) 



n = -y/Ag + 2aA2g[Re^/>(iAb/27rT) - ln(A6/27rT)] 

(B4) 
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Here, tpi^) is the digamma function. For the coherence 
(cr^(t)), one finds 
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Appendix C: Rigorous Born approximation results 
of Loss and DiVincenzo 



Since we also compare the long-time limit of our re- 
sults for the coherence {(T^{t)) with a formula derived 
in Ref. 76 by Loss and DiVincenzo using an approxima- 
tion scheme that does not apply any other approximation 
than the Born approximation. It is thus "exact" to first 
order in a. They find the steady-state value of the co- 
herence 
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/oo,LDV — 7; ^ 



A^ , Wc\/A3 



2A\ 
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(CI) 



with E = v A2 + g^j^j Euler-Mascheroni number C, 
which agrees perfectly with the Bethe ansatz prediction 
and our numerical SSE result for a < 0.1. 
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